---
title: 'Replication Codes for ''***What is Populism Good for? An Experimental Test of Mobilization Effects***'''
author: "Alexander Kustov and Yaoyao Dai"
output:
  html_notebook:
    df_print: paged
    theme: cerulean
    toc: yes
    toc_depth: 3
    toc_float:
      collapsed: no
      smooth_scroll: yes
  html_document:
    df_print: paged
    theme: cerulean
    highlight: zenburn
    toc: yes
    toc_depth: 3
    toc_float:
      collapsed: no
      smooth_scroll: yes
editor_options:
  chunk_output_type: inline
---


***
This is an R notebook with codes to replicate all the tables and figures in the paper *"What is Populism Good for? An Experimental Test of Mobilization Effects"* and its Online Appendices.To replicate the analysis, please download the replication data to the same folder as this .rmd file and run the .rmd file in R. 

# Load data and packages



```{r}
library(lmtest)
library(multiwayvcov)
library(sandwich)
library(stargazer)
library(TOSTER)

library(ggplot2)
library(estimatr)
library(dplyr)
library(tibble)
library(ggtext)
library(broom)
```

```{r}
# cleaned candidate-pair data for main analyses
load("rap_paired.RData")
# cleaned candidate level data for replication of Dai&Kustov (2022)
load("rap_conjoint_by_cand.RData")
# survey responses with original values
data <- read.csv("data_clean.csv")
```

## Codebook for key variabels in candidate-pair data

Variable  | Variable Description | Values  
-|-----|:------
voted | Whether respondent voted for either candidate in a candidate-pair |$1$: Voted for either candidate;<br> $0$: Abstained.
policy2congruent_any | Whether any candidate in a pair offers congruent policy position as the respondent, measured by expressed attitudes towards immigration and economic policies |$1$: There is one or more policy congruent candidate;<br> $0$: Neither candidate is policy congruent.
party2congruent_any | Whether any candidate in a pair offers congruent policy position as the respondent, measured by respondent's party affiliation |$1$: There is one or more congruent candidate;<br> $0$: Neither candidate is congruent.
populist_any | Whether any candidate in a pair uses populist campaign | $1$: There is one or more populist candidate;<br> $0$: Neither candidate is populist.
populismatt | Whether respondent has populist attitudes | "non-populist"; "populist"
pop_policy_conditions | Combinations of populist rhetoric and policy positions in a candidate-pair | "Both candidates non-congruent and one populist"; <br> "Both candidates non-congruent and both populist"; <br> "One candidate congruent and populist"; <br> "One candidate congruent and non-populist"; <br> "One candidate congruent and neither populist";<br> "One candidate congruent and both populist";<br> "Both candidates congruent and both non-populist";<br> "Both candidates congruent and one populist";<br> "Both candidates congruent and both populist";<br>  "Neither candidates congruent and one populist", 
***

## Codebook for key variabels in candidate level data

Variable  | Variable Description | Values  
-|-----|:------
Y | DV: whether respondent voted for the candidate |$1$: The candidate was chosen from a pair;<br> $0$: Otherwise
policy2congruent | Whether the candidate has congruent policy position as the respondent, measured by expressed attitudes towards immigration and economic policies | "non-congruent"; "congruent"
populist | Whether the candidate uses populist rhetoric | $1$: The candidate uses populist rhetoric;<br> $0$: The candidate uses liberal democratic rhetoric
populismatt | IV: whether respondent has populist attitudes | "non-populist"; "populist"

***

# Descriptive Statistics 

### Appendix A Table 1


```{r}
summary_data <- as.data.frame(data[, c("Female", "College", "Republican2", "Democrat2", "Populist2")])

stargazer(
  summary_data,
  type = "text",
  summary = TRUE,
  covariate.labels = c("Gender", "College", "Republican", "Democrat","Populist"),
  digits = 2,
  title = "Descriptive Statistics (n = 3502/28016)",
  iqr = TRUE
)
```


# Hypotheses Testing

### H0 \& H3: Table A3 \& Figure 1

```{r}
# main tests
mod_h0 <- lm(voted ~ policy2congruent_any, data = d_stack_pair_merged)
mod_h3 <- lm(voted ~ populist_any, data = d_stack_pair_merged)
# exploratory
mod_exp0 <- lm(voted ~ populist_any + policy2congruent_any, data = d_stack_pair_merged)
mod_exp1 <- lm(voted ~ populist_any + policy2congruent_any + populist_any : policy2congruent_any, data = d_stack_pair_merged)
```

```{r}
cl.se.mod_h0 <- sqrt(diag(cluster.vcov(mod_h0, d_stack_pair_merged$id)))
cl.se.mod_h3 <- sqrt(diag(cluster.vcov(mod_h3, d_stack_pair_merged$id)))
cl.se.mod_exp0 <- sqrt(diag(cluster.vcov(mod_exp0, d_stack_pair_merged$id)))
cl.se.mod_exp1 <- sqrt(diag(cluster.vcov(mod_exp1, d_stack_pair_merged$id)))
```

```{r}
# Table A3 in Appendix 
stargazer(mod_h0, mod_h3, mod_exp0, mod_exp1,
          type = "text",
          se=list(cl.se.mod_h0, cl.se.mod_h3, cl.se.mod_exp0, cl.se.mod_exp1),
          star.cutoffs = c(0.05, 0.01, 0.001),
          no.space=TRUE, column.sep.width = "-5pt", 
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          dep.var.labels = "Probability of Voting")
```

```{r}
cov_cluster_h0 <- vcovCL(mod_h0, cluster = d_stack_pair_merged$id) 
mod_h0_cluster <- coeftest(mod_h0, vcov. = cov_cluster_h0)
tidy_mod_h0 <- tidy(mod_h0_cluster, conf.int = TRUE)
```

```{r}
cov_cluster_h3 <- vcovCL(mod_h3, cluster = d_stack_pair_merged$id)
mod_h3_cluster <- coeftest(mod_h3, vcov. = cov_cluster_h3)
tidy_mod_h3 <- tidy(mod_h3_cluster, conf.int = TRUE)
```

```{r}
grab_coef <- function(tidy_mod, term) {
  est  <- tidy_mod$estimate[tidy_mod == term]
  lo <- tidy_mod$conf.low[tidy_mod == term]
  hi <- tidy_mod$conf.high[tidy_mod == term]
  list(est = est,
       lo  = lo,
       hi  = hi)
}

h0 <- grab_coef(tidy_mod_h0, "policy2congruent_any")
h3 <- grab_coef(tidy_mod_h3, "populist_any")
```

```{r}
plot_data <- tibble(
  row_id   = 1:6,
  attribute = c("Policy position (H0)", rep("", 2),
                "Populist rhetoric (H3)", rep("", 2)),
  level     = c("", "AB Non-congruent", "A OR B Congruent",
                "", "AB Non-populist",  "A OR B Populist"),
  estimate  = c(NA, 0,               h0$est,
                NA, 0,               h3$est),
  ci_lower  = c(NA, NA, h0$lo,
                NA, NA, h3$lo),
  ci_upper  = c(NA, NA, h0$hi,
                NA, NA, h3$hi)
) %>%
  mutate(
    attribute_fmt = ifelse(attribute != "",
                           paste0("<b>", attribute, "</b>"), ""),
    level_fmt     = ifelse(level != "",
                           paste0("&nbsp;&nbsp;&nbsp;", level), "")
  )
```

```{r}
x_range <- range(plot_data$ci_lower, plot_data$ci_upper, na.rm = TRUE)
pad     <- diff(x_range) * 0.05                 # small padding for labels
x_limits <- c(x_range[1] - pad, x_range[2] + pad)

ggplot(plot_data) +                                             # <- use the right data
  geom_vline(xintercept = 0, linetype = "dashed", colour = "gray50") +
  
  geom_point(aes(estimate, row_id), size = 3, na.rm = TRUE) +
  geom_errorbarh(aes(xmin = ci_lower, xmax = ci_upper, y = row_id),
                 height = 0.2, na.rm = TRUE) +
  
  ## headers: nudge down a bit
  geom_richtext(
    aes(label = attribute_fmt, x = -0.3, y = row_id),
    hjust = 0, size = 5, fill = NA, label.colour = NA,
    position = position_nudge(y = -0.30)
  ) +
  
  ## indented levels: smaller nudge
  geom_richtext(
    aes(label = level_fmt, x = -0.3, y = row_id),
    hjust = 0, size = 5, fill = NA, label.colour = NA,
    position = position_nudge(y = -0.10)
  ) +
  
  scale_y_reverse(breaks = NULL, expand = c(0.05, 0.05)) +
  scale_x_continuous(limits = c(-0.3, 0.3), breaks = seq(-0.05, 0.4, 0.05)) +
  labs(
    #title = "Effects of Candidate Pair Attributes on Probability of Voting for Either Candidate",
    #subtitle = "Estimates with 95% Confidence Intervals"
    x = "\nEstimated AMCE",
    y = "",
  ) +
  theme_minimal(base_size = 15) +
  theme(
    panel.grid.major.y = element_blank(),
    axis.text.y  = element_blank(),
    axis.ticks.y = element_blank(),
    plot.margin  = margin(0, 0, 0, 0, "pt")
  )
```

### H1 \& H2: Table A4 \& Figure 2 

```{r}
mod_h1 <- lm(voted ~ populist_congruent_any, 
             data = d_stack_pair_merged[d_stack_pair_merged$populismatt == "populist" &
                                          d_stack_pair_merged$policy2congruent_any == 1,])
mod_h2 <- lm(voted ~ populist_congruent_any, 
             data = d_stack_pair_merged[d_stack_pair_merged$populismatt == "non-populist" &
                                          d_stack_pair_merged$policy2congruent_any == 1,])
mod_int <- lm(voted ~ populist_congruent_any + populismatt + populist_congruent_any:populismatt, 
              data = d_stack_pair_merged[d_stack_pair_merged$policy2congruent_any == 1,])
```

```{r}
cl.se.mod_h1 <- sqrt(diag(cluster.vcov(mod_h1, 
                                       d_stack_pair_merged[d_stack_pair_merged$populismatt == "populist" &
                                                             d_stack_pair_merged$policy2congruent_any == 1,]$id)))
cl.se.mod_h2 <- sqrt(diag(cluster.vcov(mod_h2, 
                                       d_stack_pair_merged[d_stack_pair_merged$populismatt == "non-populist" &
                                                             d_stack_pair_merged$policy2congruent_any == 1,]$id)))
cl.se.mod_int <- sqrt(diag(cluster.vcov(mod_int, 
                                        d_stack_pair_merged[d_stack_pair_merged$policy2congruent_any == 1,]$id)))
```

```{r}
stargazer(mod_h1, mod_h2, mod_int,
          type = "text",
          se=list(cl.se.mod_h1, cl.se.mod_h2, cl.se.mod_int), 
          star.cutoffs = c(0.05, 0.01, 0.001),
          no.space=TRUE, column.sep.width = "-3pt", 
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          dep.var.labels = "Probability of Voting")
```


```{r}
## ──Grab the two marginal effects  (+/- 95 % CI) ───────────────────────────
coefs <- coef(summary(mod_int))

b_populist   <- coefs["populist_congruent_any",          "Estimate"]
se_populist  <- cl.se.mod_int["populist_congruent_any"]

b_int        <- coefs["populist_congruent_any:populismattpopulist", "Estimate"]
se_int       <- cl.se.mod_int["populist_congruent_any:populismattpopulist"]

## effect for non‑populists  = β_populist
## effect for populists      = β_populist + β_int
est_np  <- b_populist
se_np   <- se_populist

est_p   <- b_populist + b_int
se_p    <- sqrt(se_populist^2 + se_int^2 + 2 * vcov(mod_int)["populist_congruent_any",
                                                         "populist_congruent_any:populismattpopulist"])
```

```{r}
est_p <- coef(summary(mod_h1))["populist_congruent_any", "Estimate"]
est_np <- coef(summary(mod_h2))["populist_congruent_any", "Estimate"]

se_p <- cl.se.mod_h1["populist_congruent_any"]
se_np <- cl.se.mod_h2["populist_congruent_any"]
```


```{r}
## build plotting tibble
plot_dat <- tibble(
  subgroup   = c("Populist respondents (H1)", "Non-populist respondents (H2)"),
  estimate   = c(est_p, est_np),
  ci_lower   = estimate - 1.96 * c(se_p, se_np),
  ci_upper   = estimate + 1.96 * c(se_p, se_np),
  y_pos      = 1:2
)
```

```{r}
ggplot(plot_dat, aes(estimate, y_pos, colour = subgroup)) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "gray50") +
  
  ## keep legend: drop `show.legend = FALSE`
  geom_point(size = 4) +                                        # legend ON
  geom_errorbarh(aes(xmin = ci_lower, xmax = ci_upper), height = .18) +
  
  scale_y_reverse(expand = c(.3, .3)) +
  scale_x_continuous(limits = c(-0.2, 0.2),
                     breaks  = seq(-0.2, 0.2, 0.1)) +
  scale_colour_manual(
    values = c("Populist respondents (H1)"     = "black",
               "Non-populist respondents (H2)" = "gray50"),
    name   = "Respondent type"                 # legend title
  ) +
  labs(
    x = "Estimated AMCE",
    y = NULL
    #title    = "Figure 2: Effects of Using Populist Rhetoric by \n Policy Congruent Candidate on Probability of Voting for \n Either Candidate among Populist and Non-Populist Respondents (H1, H2)",
  ) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid.major.y = element_blank(),
    axis.text.y        = element_blank(),
    axis.ticks.y       = element_blank(),
    legend.position    = "bottom"              # place legend underneath
  )
```

### Exploratory: Table A5 \& Figure 3

```{r}
mod_exp_pop <- lm(voted ~ pop_policy_conditions, 
                  data = d_stack_pair_merged[d_stack_pair_merged$populismatt == "populist",])
mod_exp_npop <- lm(voted ~ pop_policy_conditions, 
                   data = d_stack_pair_merged[d_stack_pair_merged$populismatt == "non-populist",])
mod_exp_full <- lm(voted ~ pop_policy_conditions, data = d_stack_pair_merged)
```

```{r}
stargazer(mod_exp_pop, mod_exp_npop,mod_exp_full,
          type = "text",
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Populist Voter", "Non-populist Voter", "All Voters"),
          covariate.labels=c("Both candidates non-congruent and one populist", "Both candidates non-congruent and both populist", 
                             "One candidate congruent and populist", "One candidate congruent and non-populist", 
                             "One candidate congruent and neither populist",
                             "One candidate congruent and both populist", "Both candidates congruent and both non-populist", 
                             "Both candidates congruent and one populist", "Both candidates congruent and both populist",
                             "Constant"),
          no.space=TRUE, column.sep.width = "-3pt", 
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          dep.var.labels = "Probability of Voting")

```

```{r}
tidy_exp_pop <- tidy(mod_exp_pop, conf.int = TRUE)  %>% mutate(model = "Populist voters")
tidy_exp_npop  <- tidy(mod_exp_npop, conf.int = TRUE) %>% mutate(model = "Non-populist voters")
```

```{r}
coefs <- bind_rows(tidy_exp_pop, tidy_exp_npop) %>%
  filter(term != "(Intercept)") %>%                      # drop intercept
  mutate(term = factor(term, levels = rev(unique(term)))) # control order on y-axis
```

```{r}
coefs$term <- c("Neither candidates congruent and one populist", "Neither candidates congruent and both populist", 
                             "One candidate congruent and populist", "One candidate congruent and non-populist", 
                             "One candidate congruent and neither populist",
                             "One candidate congruent and both populist", "Both candidates congruent and both non-populist", 
                             "Both candidates congruent and one populist", "Both candidates congruent and both populist",
                 "Neither candidates congruent and one populist", "Neither candidates congruent and both populist", 
                             "One candidate congruent and populist", "One candidate congruent and non-populist", 
                             "One candidate congruent and neither populist",
                             "One candidate congruent and both populist", "Both candidates congruent and both non-populist", 
                             "Both candidates congruent and one populist", "Both candidates congruent and both populist"
                 )
```

```{r}
coefs <- coefs[, c("term","estimate", "conf.low", "conf.high", "model")]
baseline <- data.frame(
  term = c("Neither candidate congruent and neither populist", "Neither candidate congruent and neither populist"), estimate = c(0, 0),
  conf.low = c(0, 0), conf.high = c(0, 0), model = c("Populist voters", "Non-populist voters")
)
coefs <- rbind(baseline, coefs)

```

```{r}

ggplot(coefs, aes(x = estimate, y = term, color = model)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_point(position = position_dodge(width = 0.5), size = 2.6) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                 position = position_dodge(width = 0.5), height = 0) +
  scale_color_manual(values = c("Populist voters" = "black", "Non-populist voters" = "gray50")) +
  labs(x = "Estimated AMCE", y = NULL, color = NULL) +
  theme_minimal(base_size = 12)
```


# Robustness Checks 

### H0 \& H3: Table A6

```{r}
#robustness check party-congruence
mod_h0_party <- lm(voted ~ party2congruent_any, data = d_stack_pair_merged)
mod_h3_party <- lm(voted ~ populist_any, data = d_stack_pair_merged)
mod_exp0_party <- lm(voted ~ populist_any + party2congruent_any, data = d_stack_pair_merged)
mod_exp1_party <- lm(voted ~ populist_any + party2congruent_any + populist_any : party2congruent_any, data = d_stack_pair_merged)
```

```{r}
cl.se.mod_h0_party <- sqrt(diag(cluster.vcov(mod_h0_party, d_stack_pair_merged$id)))
cl.se.mod_h3_party <- sqrt(diag(cluster.vcov(mod_h3_party, d_stack_pair_merged$id)))
cl.se.mod_exp0_party <- sqrt(diag(cluster.vcov(mod_exp0_party, d_stack_pair_merged$id)))
cl.se.mod_exp1_party <- sqrt(diag(cluster.vcov(mod_exp1_party, d_stack_pair_merged$id)))
```

```{r}
stargazer(mod_h0_party, mod_h3_party, mod_exp0_party, mod_exp1_party,
          type = "text",
          se=list(cl.se.mod_h0_party, cl.se.mod_h3_party, cl.se.mod_exp0_party, cl.se.mod_exp1_party),
          star.cutoffs = c(0.05, 0.01, 0.001),
          no.space=TRUE, column.sep.width = "-5pt", 
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          dep.var.labels = "Probability of Voting")
          
```


### H1 \& H2: Table A7

```{r}
# robustness check party-congruence 
mod_h1_party <- lm(voted ~ populist_partycongruent_any, data = d_stack_pair_merged[d_stack_pair_merged$populismatt == "populist" & d_stack_pair_merged$party2congruent_any == 1,])
mod_h2_party <- lm(voted ~ populist_partycongruent_any, data = d_stack_pair_merged[d_stack_pair_merged$populismatt == "non-populist" & d_stack_pair_merged$party2congruent_any == 1,])
mod_int_party <- lm(voted ~ populist_partycongruent_any + populismatt + populist_partycongruent_any:populismatt, data = d_stack_pair_merged[d_stack_pair_merged$party2congruent_any == 1,])
```

```{r}
stargazer(mod_h1_party, mod_h2_party, mod_int_party,
          type = "text",
          star.cutoffs = c(0.05, 0.01, 0.001),
          no.space=TRUE, column.sep.width = "-3pt", 
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          dep.var.labels = "Probability of Voting")
```

# Replication of Dai and Kustov (2022): Table A2

```{r}
mod_DK_pop <- lm(Y ~ populist + policy2congruent, data = d_stack_cand_merged[d_stack_cand_merged$populismatt == "populist",])
mod_DK_nonpop <- lm(Y ~ populist + policy2congruent, data = d_stack_cand_merged[d_stack_cand_merged$populismatt == "non-populist",])
mod_DK_int <- lm(Y ~ populist + policy2congruent + populismatt + populist:populismatt, data = d_stack_cand_merged)
```

```{r}
cl.se.mod_DK_pop <- sqrt(diag(cluster.vcov(mod_DK_pop, d_stack_cand_merged[d_stack_cand_merged$populismatt == "populist",]$id)))
cl.se.mod_DK_nonpop <- sqrt(diag(cluster.vcov(mod_DK_nonpop, d_stack_cand_merged[d_stack_cand_merged$populismatt == "non-populist",]$id)))
cl.se.mod_DK_int <- sqrt(diag(cluster.vcov(mod_DK_int, d_stack_cand_merged$id)))
```

```{r}
stargazer(mod_DK_pop, mod_DK_nonpop, mod_DK_int,
          type = "text",
          se=list(cl.se.mod_DK_pop, cl.se.mod_DK_nonpop, cl.se.mod_DK_int), 
          star.cutoffs = c(0.05, 0.01, 0.001),
          no.space=TRUE, column.sep.width = "-10pt", 
          omit.stat=c("rsq", "ser", "f"), notes = "", notes.append = F,  notes.label = "",
          dep.var.labels = "Probability of Selection")

```


